package fem2.material.cohesive_law2;

import inf.jlinalg.Array2DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;
import java.util.Map;
import java.util.TreeMap;

import math2.MathUtilities;
import math2.MutableDouble;
import fem2.Element;
import fem2.Variable;
import fem2.element.crack_sda2.SDAElement;
import fem2.enu.VariableName;
import fem2.jlinalg.DiagonalMatrix;
import fem2.material.FluencyCriteria;

/**
 * Mode I cohesive law based on Sancho paper: An embedded crack model for finite
 * element analysis of concrete fracture.<br>
 * Note: this class works with constant opening only
 * 
 * @author hbui
 * 
 */
public class SanchoCohesiveLaw extends CohesiveLaw {

	private Map<Element, Boolean> isLoading = new TreeMap<Element, Boolean>();

	private FluencyCriteria softeningRule;

	public SanchoCohesiveLaw(double Gf, double ft, FluencyCriteria softeningRule) {
		super(Gf, ft);
		this.softeningRule = softeningRule;
	}

	/**
	 * compute value of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected double computeSofteningLaw(double x) {
		return softeningRule.computeFluencyValue(x);
	}

	/**
	 * compute value of derivative of softening law
	 * 
	 * @param x
	 *            opening norm
	 * @return
	 */
	protected double computeSofteningLawDerivative(double x) {
		return softeningRule.computeFluencyDerivative(x);
	}

	@Override
	public IVector getJumpTraction(Element e, int k) {
		int dim = e.getMeshPart().getDimension();
		IVector w = getVariable(e).get(k).getVectorVariable(VariableName.DISCONTINUITY_JUMP);
		IVector t = (IVector) w.clone();

		/*
		 * check for zero opening
		 */
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);
		if (wtilde < MathUtilities.tol6) {
			double ft = computeSofteningLaw(0);
			t = new ArrayVector(dim);
			t.set(0, ft);
			return t;
		}

		double wt_old = getVariable(e).get(k).getScalarVariable(VariableName.EQUIVALENT_JUMP)
				.getValue();

		if (!isLoading.get(e)) {

			/*
			 * unloading
			 */
			double coeff = computeSofteningLaw(wt_old) / wt_old;
			BLAM.scale(coeff, t);

		} else {

			/*
			 * loading
			 */
			double coeff = computeSofteningLaw(wtilde) / wtilde;
			BLAM.scale(coeff, t);

		}

		return t;
	}

	@Override
	public IMatrix getJumpStiffness(Element e, int k) {
		int dim = e.getMeshPart().getDimension();
		IVector w = getVariable(e).get(k).getVectorVariable(VariableName.DISCONTINUITY_JUMP);
		IMatrix C;

		/*
		 * check for zero opening
		 */
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);
		if (wtilde < MathUtilities.tol6) {

			C = new DiagonalMatrix(dim);
			// BLAM.zero(C);

			// C = new Array2DMatrix(dim, dim);
			// C.set(0, 0, 1e10);
			// C.set(0, 1, -1e10);
			// C.set(1, 0, -1e10);
			// C.set(1, 1, 1e10);

			return C;
		}

		double wt_old = getVariable(e).get(k).getScalarVariable(VariableName.EQUIVALENT_JUMP)
				.getValue();

		if (!isLoading.get(e)) {
			/*
			 * unloading
			 */
			C = new DiagonalMatrix(dim);
			double coeff = computeSofteningLaw(wt_old) / wt_old;
			DiagonalMatrix.set(C, coeff);

			// Debugger.watch(e.getId() + ", wt_old (t) = " + wt_old);

		} else {

			/*
			 * loading
			 */
			C = new Array2DMatrix(dim, dim);
			double f = computeSofteningLaw(wtilde);
			double df = computeSofteningLawDerivative(wtilde);
			if (dim == 2) {
				double wn = w.get(0);
				double wm = w.get(1);
				double tmp1 = f / wtilde;
				double tmp2 = Math.pow(wn / wtilde, 2);
				double tmp3 = Math.pow(wm / wtilde, 2);
				double tmp4 = wn * wm / (wtilde * wtilde);
				double c00 = tmp1 * tmp3 + df * tmp2;
				double c01 = -tmp1 * tmp4 + df * tmp4;
				double c11 = tmp1 * tmp2 + df * tmp3;
				C.set(0, 0, c00);
				C.set(0, 1, c01);
				C.set(1, 0, c01);
				C.set(1, 1, c11);
			} else if (dim == 3) {
				double wn = w.get(0);
				double wm1 = w.get(1);
				double wm2 = w.get(2);
				double tmp1 = f / wtilde;
				double tmp2 = Math.pow(wn / wtilde, 2);
				double tmp3 = Math.pow(wm1 / wtilde, 2);
				double tmp4 = Math.pow(wm2 / wtilde, 2);
				double tmp5 = wn * wm1 / (wtilde * wtilde);
				double tmp6 = wn * wm2 / (wtilde * wtilde);
				double tmp7 = wm1 * wm2 / (wtilde * wtilde);
				double c00 = tmp1 * (tmp3 + tmp4) + df * tmp2;
				double c01 = -tmp1 * tmp5 + df * tmp5;
				double c02 = -tmp1 * tmp6 + df * tmp6;
				double c11 = tmp1 * (tmp2 + tmp4) + df * tmp3;
				double c12 = -tmp1 * tmp7 + df * tmp7;
				double c22 = tmp1 * (tmp2 + tmp3) + df * tmp4;
				C.set(0, 0, c00);
				C.set(0, 1, c01);
				C.set(0, 2, c02);
				C.set(1, 0, c01);
				C.set(1, 1, c11);
				C.set(1, 2, c12);
				C.set(2, 0, c02);
				C.set(2, 1, c12);
				C.set(2, 2, c22);
			}

			// Debugger.watch(e.getId() + ", wtilde (C) = " + wtilde);
		}

		return C;
	}

	@Override
	public void Initialize(Element e) {
		int ngc = ((SDAElement) e).getNumGaussPointOnCrack();
		int dim = e.getMeshPart().getDimension();

		ArrayList<Variable> vars = new ArrayList<Variable>();

		for (int i = 0; i < ngc; i++) {
			Variable v = new Variable(0);

			v.registerVectorVariable(VariableName.DISCONTINUITY_JUMP, dim, false);
			v.registerScalarVariable(VariableName.EQUIVALENT_JUMP, false);

			vars.add(v);
		}

		setVariable(e, vars);
		isLoading.put(e, false);
	}

	@Override
	public void InitializeSolutionStep(Element e) {
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		int ngc = ((SDAElement) e).getNumGaussPointOnCrack();
		IVector w = ((SDAElement) e).getCrackOpeningParameter();
		double wtilde = BLAM.norm(BLAM.NORM_TWO, w);

		for (int k = 0; k < ngc; k++) {

			/*
			 * save the discontinuity jump
			 */
			BLAM.copy(w, getVariable(e).get(k).getVectorVariable(VariableName.DISCONTINUITY_JUMP));

			/*
			 * check for loading/unloading
			 */
			MutableDouble wt = getVariable(e).get(k)
					.getScalarVariable(VariableName.EQUIVALENT_JUMP);

			if (wtilde > wt.getValue()) {
				wt.setValue(wtilde);
				isLoading.put(e, true);
			} else {
				isLoading.put(e, false);
			}
		}

	}

	@Override
	public void FinalizeSolutionStep(Element e) {
		for (Variable v : getVariable(e)) {
			v.save();
		}
	}

	@Override
	public void CalculateOutputData(Element e) {
	}

}
